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Abstract 

The late-time entropy production by the massive particle decay induces the 
various cosmological effects in the early epoch and modify the standard sce- 
nario. We investigate the thermalization process of the neutrinos after the 
entropy production by solving the Boltzmann equations numerically. We 
find that if the large entropy are produced at t ~ 1 sec, the neutrinos are 
not thermalized very well and do not have the perfect Fermi-Dirac distri- 
bution. Then the freeze-out value of the neutron to proton ratio is altered 
considerably and the produced light elements, especially 4 He, are drastically 
changed. Comparing with the observational light element abundances, we 
find that Tr < 0.7 MeV is excluded at 95 % C.L. We also study the case in 
which the massive particle has a decay mode into hadrons. Then we find that 
Tr should be a little higher, i.e. Tr > 2.5 MeV - 4 MeV, for the hadronic 
branching ratio = 10~ 2 — 1. Possible influence of late-time entropy pro- 
duction on the large scale structure formation and temperature anisotropics 
of cosmic microwave background is studied. It is expected that the future 
satellite experiments (MAP and PLANCK) to measure anisotropies of cosmic 
microwave background radiation temperature can detect the vestige of the 
late-time entropy production as a modification of the effective number of the 
neutrino species N° s . 

98.80.Cq, 98.70.Vc, KUNS-1639 
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I. INTRODUCTION 



In the standard big bang cosmology it had been assumed tacitly that the universe was 
dominated by the thermal radiation at the early epoch. Even in the paradigm of the modern 
cosmology it is commonly believed that thermal radiation was produced by the reheating 
process after the primordial inflation and they dominated the energy of the universe at 
sufficiently early epoch. Here we ask, "How early should the universe be dominated by 
radiation in order to success the standard big bang cosmology?". We could say that the 
energy of the universe should be dominated by the radiation at least before the beginning 
of the big bang nucleosynthesis (BBN) epoch. In this paper we answer the above question. 

The various models of the modern particle physics beyond the standard model predicts 
a number of unstable massive particles which have long lifetimes and decays at about BBN 
epoch. The energy density of the non-relativistic particles or the oscillation energy density of 
the scalar fields (inflaton and so on) decreases as PNii(t) oc a(t)~ 3 , where a(t) is a scale factor. 
On the other hand since the radiation energy density decreases more rapidly p(t) oc a(t)~ A , 
if the energy density of the massive non-relativistic particles or the oscillating scalar fields 
is large enough, it immediately dominates the universe as it expands, and the universe 
necessarily becomes matter-dominated until the cosmic time reaches to their lifetime. When 
the particles decay into standard particles (e.g. photon and electron), they produce the 
large entropy and the universe becomes radiation-dominated again. It is expected that such 
process would change the initial condition for the standard big bang scenario. We call the 
process "late-time entropy production" . 

Now we have some interesting candidates for late-time entropy production in models 
based on supersymmetry (SUSY). It is known that gravitino and Polonyi field which exist 
in local SUSY (i.e. supergravity ) theories have masses of ~ (9(100GeV — lOTeV) 0. In 
addition they have long lifetimes because they interact with the other particle only through 
gravity. For example, since Polonyi field @ which has a heavy mass ~ 10 TeV cannot be 
diluted by the usual inflation, it immediately dominates the universe and decays at the BBN 
epoch. Moreover it is also known that in the superstring theories there exist many light fields 
called dilaton and moduli which have similar properties to the Polonyi field. 

Recently Lyth and Stewart || considered a mini-inflation called "thermal inflation" 
which dilutes the above dangerous scalar fields. In the thermal inflation scenario, however, 
the flaton field which is responsible for the thermal inflation decays at late times. In par- 
ticular, if Polonyi (moduli) mass is less than ~ 1 GeV which is predicted in the framework 
of gauge-mediated SUSY breaking models [[| , the sufficient dilution requires that the flaton 
decays just before BBN ||. Thus, in thermal inflation models, one should take care of the 
late-time entropy production. 

To keep the success of BBN, any long-lived massive particles or the coherent oscillation 
of any scalar fields which dominate the universe at that time must finally decay into the 
standard particles before the beginning of BBN. Moreover the decay products would have 
to be quickly thermalized through scatterings, annihilations, pair creations and further de- 
cays and make the thermal bath of photon, electron and neutrinos. Concerning photons 
and electrons which electromagnetically interact, the interaction rate is much more rapid 
than the Hubble expansion rate at that time. Therefore it is expected that the photon 
and electron which are produced in the decay and subsequent thermalization processes are 
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efficiently thermalized. The problem is that neutrinos can interact only through the weak 
interaction. In the standard big bang cosmology the neutrinos usually decouple from the 
electromagnetic thermal bath at about T ~ 2 — 3MeV. Therefore it is approximately inferred 
that they can not be sufficiently thermalized at the temperature T < a few MeV. Namely 
the reheating temperature after the entropy production process should be high enough to 
thermalize the neutrinos. Though people had ever used the rough constraints on the reheat- 
ing temperature between lMeV - lOMeV, in the previous paper |3j we pointed out that the 
neutrino thermalization is the most crucial for the successful BBN. In this paper we describe 
the detail of the method to obtain the neutrino spectrum and the formulations to integrate 
a set of Boltzmann equations numerically , and we study the constraint on the reheating 
temperature using the obtained neutrino spectrum and the full BBN network calculations 
with the revised observational light element abundances. 

The above constraint is almost model-independent and hence conservative because we 
only assume that the massive particle decay produces the entropy. However, a more stringent 
constraint can be obtained if we assume a decay mode into quarks or gluons. In this case some 
modifications are needed for the above description. When the high energy quark- ant iquark 
pairs or gluons are emitted, they immediately fragment into a lot of hadrons (pions , kaons, 
protons, neutrons, etc.). It is expected that they significantly influence the freeze-out value 
of neutron to proton ratio at the beginning of BBN through the strong interaction with the 
ambient protons and neutrons. In the previous paper || we did not consider such hadron 
injection effects on BBN. Therefore we carefully treat the hadron injection effects in the 
present paper. 

For another constraint, the late-time entropy production may induce the significant ef- 
fects on the anisotropies of the cosmic microwave background radiations (CMB). Lopez et 
al. [0] pointed out that the CMB anisotropies are very sensitive to the equal time of matter 
and radiation. When the reheating temperature is so low that neutrinos do not be suffi- 
ciently thermalized, the radiation density which consists of photon and neutrinos becomes 
less than that in the standard big bang scenario. It may give distinguishable signals in the 
CMB anisotropies as a modification of the effective number of neutrino species N„. With 
the present angular resolutions and sensitivities of COBE observation [^3| it is impossible to 
set a constraint on but it is expected that future satellite experiments such as MAP 
and PLANCK |10[ will gives us a useful information about Nff. In addition the above effect 
may also induce the signals in the observed power spectrum of the density fluctuation for 
the large scale structure as a modification of the epoch of the matter-radiation equality. 

The paper is organized as follows. In Sec. II we introduce the formulation of the basic 
equations and the physical parameters. In Sec. Ill we briefly review the current status of 
the observational light element abundances. In Sec. IV we study the spectra of the electron 
neutrino and the mu(tau)-neutrino by numerically solving the Boltzmann equations, and 
the constraints from BBN are obtained there. In Sec V we investigate the additional effects 
in the hadron injection by the massive particle decay. In Sec. VI we consider the another 
constraints which come from observations for large scale structures and anisotropies of CMB. 
Sec VII is devoted to conclusions. In Appendix we introduce the method of the reduction 
for the nine dimension integrals into one dimension. 
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II. FORMULATION OF NEUTRINO THERM ALIZATION 



A. Reheating Temperature 

In order to discuss the late-time entropy production process, we should formulate the 
equations which describe the physical process. Here the reheating temperature Tr is an 
appropriate parameter to characterize the late-time entropy production. We define the 
reheating temperature Tr by 

F = 3H(T R ), (1) 

where T is the decay rate (=r _1 ) and H(Tr) is the Hubble parameter at the decay epoch (t 
= r). Q The Hubble parameter is expressed by 

V 90 j Mq k ' 

where is the statistical degrees of freedom for the massless particles and Mq is the reduced 
Plank mass (= 2.4 x 10 18 GeV). Then the reheating temperature is given by 



T R = 0.554 v /rM G . (3) 

Here we have used g* = 43/4. From Eq. (Q), we can see that the reheating temperature has 
the one to one correspondence with the lifetime of the parent massive particle. 

Here we define the effective number of neutrino species N„ as a parameter which char- 
acterize the time evolution of the energy density of neutrinos. Here N^ s is defined by 



Pstd 

where p st d is the total neutrino energy density in the standard big bang model (i.e. no 
late-time entropy production and three neutrino species). 



B. Basic Equations 

When a massive particle which is responsible for the late-time entropy production de- 
cays, all emitted particles except neutrinos are quickly thermalized and make a thermal bath 
with temperature ~ Tr. For relatively low reheating temperature Tr < lOMeV neutrinos 
are slowly thermalized. Since in realistic situations the decay branching ratio into neutrinos 
is very small, we assume that neutrinos are produced only through annihilations of electrons 



Since the actual decay is not instantaneous, the matter-dominated universe smoothly changes into 
radiation-dominated one. Thus it is rather difficult to clearly identify the reheating temperature by 
observing the evolution of the cosmic temperature. Instead we "define" the reheating temperature 
formally by Eq. (|l]) 
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and positrons, i.e. e + + e — > Vi + Ui(i = e, /i, r). The evolution of the distribution function 



fi of the neutrino V{ is described by the momentum dependent Boltzmann equation | 1 



where the right hand side is the total collision term. Q When the reaction is two bodies 
scattering 1 + 2 — > 3 + 4, it is given by the expression, 



Ci,coll — „ j-, 



d 3 p 2 d 3 pz d 3 p<± 



2E X ^ J 2E 2 (2tt) 3 2E 3 (2it) 3 2£ 4 (2tt) 3 

x(2w) 4 6^(p 1 +p 2 -p 3 -p 4 )A(f 1 ,f 2 J 3 J,)S\M\ 2 12 _, 34 , (6) 

where \M\ 2 is the scattering amplitude summed over spins of all particles, S is the sym- 
metrization factor which is 1/2 for identical particles in initial and final states, A = 
/ 3 / 4 (l — /i)(l — f 2 ) — ^1/2(1 — /3XI — fi) is the phase space factor including Pauli blocking 
of the final states. Then the total collision term C ijCO \\ is expressed by, 

Ci,coll Cj,ann Cj,scat) (7) 

where Cj ann is the collision term for annihilation processes and Cj scat is collision term for 
elastic scattering processes. Here we consider the following processes: 

ft + ^ <-> e + + e", 
Vi + e <-> z/j + e . 

In this paper we have treated neutrinos as Majorana ones (i.e., 1/ — u). It should be noted 
that there are no differences between Majorana neutrinos and Dirac ones as long as they 
are massless, and since the temperature is 0( MeV) at least in this situation, we could have 
treated them as if they were massless particles. The relevant reactions are presented in 
Table | for v e and Table |I| for and v r . f\ 



The integrated Boltzmann equation [12| is not adequate in the present problem. As we show in 
Sec. [V, the spectral shape of the momentum distribution obtained by our scheme is much different 
from the equilibrium one. It should be noticed that the integrated Boltzmann equation assumes 
that the shape of the momentum distribution is the same as the equilibrium one. Thus we should 
solve the momentum dependent Boltzmann equation. 

3 Here we neglect the neutrino self- inter actions. It may lead to underestimate the kinetic equi- 
librium rate for high reheating temperatures. However, we think that this effect does not change 
the results very much. The interactions between electrons and neutrinos are the most important 
because they transfer the energy of the thermal bath to neutrinos. The self-interactions of the 
neutrinos cannot increase the energy density of neutrinos but mainly change their momentum dis- 
tribution. Furthermore, the neutrino number densities are much smaller than the electron number 
density at low reheating temperature with which we are concerned. Thus differences caused by the 
neutrino self interactions are expected to be small. 
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The collision terms are quite complicated and expressed by nine dimensional integrations 
over momentum space. However, if we neglect electron mass and assume that electrons 
obey the Boltzmann distribution e~ p ^ T , the collision terms are simplified to one dimensional 
integration form. f] Then C i ann is given by |T5|Jl6[ 



a 



j,ann 



2tt 2 



Pi dp'ii^iifMfiip'i) - feq{Pi)feq{p'i))i 



(8) 



where f eq {= l/(e Pi//T + 1)) is the equilibrium distribution and (av)i is the differential cross 
sections given by 



lav) 



4G| 
9tt 



(C 2 V + C\)pp\ 



AGl , ~ 2 ~ 2 



[(TV 



fl,T 



9vr 



(Cy + C/)pp', 



(9) 

(10) 
(11) 



where we take Cy — \ + 2 sin 2 0w, Ca = |, Cy = Cy — 1 {Ca = Ca — 1) and 0w is Weinberg 
angle (sin 2 9 W ~ 0.231) 0. 

As for elastic scattering processes, C itScat is also simplified to one dimensional integration 
(see Appendix), and it is expressed as 



a 



G 



i,scat 



2tt 



F (Cy + CD 



"4 ( L MFiiPhPdi 1 - hip'i)) + I dpiFziPhP'di 1 - hip'i) 



Pi VJo 



1 ~ fiiPi 
Pi 



(12) 



~ 2 ~ 2 

where (C v + C^) should be replaced by (Cy + Ca ) for i = /x, r, and the functions 
F lj F 2 ,B ll B 2 are given by 



F 1 (p,p') = D(p,p') + E(p,p')e- p '/ T , 
F 2 {p,p') = D(p',p)e^- p ' )/T + E(p,p')e- p '/ T , 
Bi(p,p) = F 2 (p',p), B 2 (p,p') = Fi(p',p), 



(13) 



where 



4 The errors due to neglecting the electron mass is small and the deviation is just a few percent. 
We show the reasons as follows. The difference between Fermi-Dirac and Maxwell-Boltzmann 
distribution "df" is less than one at most df < 1.0. The week interaction rate is almost expressed 
by (av)n e /H(t), where (av) ~ G\vr? e and n e is an electron number density. Then the error is at 
most estimated by, (awv}n e /H(t) x df ^ 10~ 2 (for T < 0.5MeV). Therefore the deviation is a few 
percent and the neglecting the electron mass does not change the results. The other methods of 
the approximation to reduce the integral from nine to two dimensions in which the electron mass 
is not neglected are presented in ref. [13,14] 



6 



D(p,p') = 2T\p 2 + p' 2 + 2T{p - p') + AT 2 ), 
E(p,p') = -T 2 [p 2 p' 2 + 2pp'(p + p')T 

+2(p + p') 2 T 2 + 4(p + p')T 3 + 8T 4 ]. (14) 

Together with the above Boltzmann equations, we should solve the energy-momentum 
conservation equation in the expanding universe: 

^j± = -3H(t)(p(t) + P(t)), (15) 

where p(t) = p ( f ) + p 7 + p e + p u is the total energy density of 0, photon, electron and neutrinos 
and it is given by 

, (t ) = H(t) + =3 + + ^/<%V,( 5 ) + (16) 



where _E e = yg 2 + m 2 is the electron energy. P(t) = P 7 (t) + P e ±(t) + P v {t) is the total 
pressure, 

P « = ^ + l/ 3g .[e,p(g rT ) + 1 ] + ^/^fa> + J^/ < 17 > 
(t) is the Hubble parameter, 

H{() = W) = Tsk^- (18) 

The time evolution equation of p^ is given by 

^± = -r P4> -3H Pti> . (19) 

Practically we solve the time evolution of the photon temperature instead of Eq. ([15]), 

cff 7 _ -p^/r^ + AHp 1 + 3H(p e ± + P e ±) + 4ifp^ + dp u /dt 



dt dp-y I <9T 7 \ a (t)+ dp e ± 1 <9T 7 1 a(t) 

together with Eqs. (§), (JT|) and (pp. 



(20) 



III. OBSERVATIONAL LIGHT ELEMENT ABUNDANCES 

In this section we briefly show the current status of the observational light element 
abundances. Concerning the deuterium abundance, the primordial D/H is measured in the 
high redshift QSO absorption systems. For the most reliable D abundance, we adopt the 
following value which is obtained by the clouds at z = 3.572 towards Q1937-1009 and at z 
= 2.504 towards Q1009+2956 0, 

D/H = (3.39 ±0.25) x 10" 5 . (21) 
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On the other hand, recently the high deuterium abundance is reported in relatively low 
redshift absorption systems at z = 0.701 towards Q1718+4807 @, D/H = (2.0±0.5) x 10~ 4 . 
Another group also observes the clouds independently JOJ. However, because they do not 
have full spectra of the Lyman series, the analyses would be unreliable. More recently 



Kirkman et al. |20| observed the quasar absorption systems at z = 2.8 towards Q0130-4021 
and they obtain the upper bound, D/H < 6.7 x 10 -5 . Moreover Molaro et al. reported D/H 
~ 1.5 x 10~ 5 which was observed in the absorber at z = 3.514 towards APM 08279+5255 



although it has the large systematic errors in the hydrogen column density |2T| . Considering 



the current situation, we do not adopt the high deuterium value in this paper. 

The primordial 4 He mass fraction Y p is observed in the low metalicity extragalactic HII 
regions. Since 4 He is produced with Oxygen in the star, the primordial value is obtained 
to regress to the zero metalicity O/H — > for the observational data. Using the 62 blue 
compact galaxies (BCG) observations, it was reported that the primordial Y is relatively 
" low", Y p ~ 0.234 |[22|| . However, recently it is claimed that Hel stellar absorption is an 
important effect though it was not included in the previous analysis [^] properly. They 
found the relatively "high" primordial value, Y p = 0.245 ± 0.004. More recently Fields and 
Olive |M| also reanalyze the data including the Hel absorption effect and they obtain 

Y p = 0.238 ± (0.002) stat ± (0.005)^, (22) 

where the first error is the statistical uncertainty and the second error is the systematic one. 
We adopt the above value as the observational Y p . 

The primordial 7 Li/H is observed in the Pop II old halo stars. In general a halo star 
whose surface effective temperature is low (the mass is small), has the deep convective zone. 
For such a low temperature star, the primordial 7 Li is considerably depleted in the warm 
interior of the star. On the other hand for the high temperature stars (T e ff > 5500K), it is 
known that the primordial abundance is not changed and they have a "plateau" of the 7 Li 
as a function of the effective temperature. In addition, though it is also known that 7 Li/H 
decreases with decreasing Fe/H, 7 Li still levels off at lower metalicity, [Fe/H]< —1.5, in 
the plateau stars. We adopt the recent measurements which are observed by Bonifacio and 



Molaro |25[ . They observed 41 old halo stars which have the plateau. We take the additional 
larger systematic error, because there may be underestimates in the stellar depletion and 
the production by the cosmic ray spallation. Then we obtain 

log 10 ( 7 Li/H) = -9.76 ± (0.012) stat ± (0.05) syst ± (0.3) add . (23) 



IV. NEUTRINO THERM ALIZATION AND BBN 

A. Time evolution of Neutrino spectrum 

The evolution of the cosmic temperature T is shown in Fig. [I] (a) for Tr = 10 MeV, 
and (b) for Tr = 2 MeV. In Fig. [I] (a), it is seen that the temperature decreases slowly 
as t~ 1//4 , i.e. a -3 / 8 before the decay epoch, t ~ T _1 (~ 5 x 10~ 2 sec) which corresponds to 
T R = 10 MeV. This is because the actual decay is not instantaneous and <p decays into 
radiation continuously at the rate V ||26|| . Then The universe is still in M.D. After the decay 
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epoch t ^> r _1 , all 0-particles decay and the temperature decreases as a -1 and t~ l l 2 . Then 
the universe becomes radiation-dominated epoch. Since at the temperature T < 0.5MeV 
(t ^ 3 sec), electrons and positrons annihilate into photons e + e _ — > 27, the temperature is 
slightly heated. From Fig. [T| (b) we can see that the temperature decreases as t _1//4 until 
the decay epoch (t < 0.1 sec) which corresponds to Tr = 2 MeV. After the decay epoch, 
the temperature decreases as i -1 / 2 (R.D.). In the actual computation we take a initial 
condition that there exists the net radiation energy density though the universe is in M.D. 
This represents the situation that the massive particle necessarily dominated the universe as 
it expands. On the other hand even if there are at first no radiation p R ~ 0, i.e.T ~ which 
corresponds to the initial condition of the oscillation epoch after the primordial inflation or 
thermal inflation, the cosmic temperature immediately tracks the same curve t -1 / 4 and then 
their decay establish the radiation dominated universe T oc t -1 / 2 . Therefore our treatment 
is quite a general picture for each entropy production scenario and it does not depend on 
whether the net initial radiation energy exists or not, only if once the unstable non-relativistic 
particles dominate the energy density of the universe. 

In Fig. ^|we show the evolutions of p Ue and p u (=/v) (a) for T R = 10 MeV and (b)2 MeV. 
From Fig. |2|(a) we can see that if Tr = 10 MeV, cosmic energy density is as same as the 
case of standard big bang cosmology. As shown in Fig. 0(b), however, the energy density 
of each neutrino species for Tr = 2 MeV is smaller than the case of standard scenario. 
Since the electron neutrinos interact with electrons and positrons through both charged and 
neutral currents, they are more effectively produced from the thermal bath than the other 
neutrinos which have only neutral current interactions. The final distribution functions f e 
and / M (= f T ) are shown in Fig. g (a) for T R = 10 MeV and (b) 2 MeV. For T R = 10 MeV, 
each neutrino is thermalized well and the perfect Fermi-Dirac distribution is established. 
For Tr = 2 MeV, however, the distributions are not thermal equilibrium forms. [] 



5 As we noted in Sec. [Oj, we must not use the integrated Boltzmann equation instead of the 
momentum dependent Boltzmann equation in the present problem because the former assumes 
the equilibrium distribution. To see this, let us define the ratio Re for a neutrino species by 
Re = (Pu/n u )/(3.151T u ), where p u is the neutrino energy density, n v is the neutrino number 
density, T u is the effective neutrino temperature which is defined by the neutrino number density 

~ 1/3 

as, T v = (27r 2 /(3C(3))n^) . Here both p v and n v are computed by integrating the neutrino 
distribution function which is obtained by solving the momentum dependent Boltzmann equation. 
Approximately Re represents a ratio of the mean energy per a neutrino to the thermal equilibrium 
one. If the neutrino is in thermal equilibrium, Re is unity. In the case of the integrated Boltzmann 
equation, because it is assumed that the shape of the neutrino distribution is the same as the 
equilibrium one at any times, Re is necessarily unity. On the other hand, in the case of our 
scheme, i.e. the momentum dependent Boltzmann equation, Re can not be unity. We have 
computed the ratio Re in some representative reheating temperatures for electron neutrino and 
have found that they deviated from unity more at the lower reheating temperatures, Re = 1.00, 
1.03 and 1.50 (for Tr = 10 MeV, 3 MeV and 1 MeV). Moreover at the lower reheating temperature 
than 1 MeV, the deviation is much larger. This result tells us that the neutrino distribution 
deviates from the thermal equilibrium shape considerably at the low reheating temperatures and 
we should solve the momentum dependent Boltzmann equation. Re has a tendency to increase as 
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In Fig. |] we can see the change of the effective number of neutrino species Nff as a 
function of the reheating temperature Tr. If Tr > 7 MeV, N„ is almost equal to three 
and neutrinos are thermalized very well. We can regard that it corresponds to the initial 
condition which has ever been used for the standard big bang cosmology. On the other 
hand, if Tr ^$ 7 MeV, becomes smaller than three. 



B. Neutrino thermalization and neutron to proton ratio 

If the neutrinos are not thermalized sufficiently and do not have the perfect Fermi-Dirac 
distribution, i. e. in this case there is the deficit of the neutrino distribution due to the low 
reheating temperature, it considerably influences the produced light element abundances. 
In particular, the abundance of the primordial 4 He is drastically changed. The change 
of the neutrino distribution function influences the neutrino energy density and the weak 
interaction rates between proton and neutron. At the beginning of BBN (T ~ 1 MeV - 0.1 
MeV) the competition between the Hubble expansion rate H and the weak interaction rates 
T n ^ p determines the freeze-out value of neutron to proton ratio n/p. After the freeze-out 
time, neutrons can change into protons only through the free decay with the life time r n . 
Since 4 He is the most stable light element and the almost all neutrons are synthesized into 
4 He, the abundance of the primordial 4 He is sensitive to the freeze-out value of neutron to 
proton ratio. 

If the neutrino energy density gets smaller than that of the standard BBN (SBBN), 
Hubble expansion rate which is proportional to the square of the total energy density is 
also decreased. Then the freeze out time becomes later and the (3 equilibrium between 
neutrons and protons continues for longer time. As a result less neutrons are left. In this 
case the predicted 4 He is less than the prediction of SBBN. The effect due to the speed-down 
expansion is approximately estimated by 

Ay ~ -0.1{-A Ptot / Ptot ), (24) 

where Y is the mass fraction of 4 He and ptot is the total energy density of the universe. 

Moreover, when the electron neutrino is not thermalized, there is interesting effect by 
which more 4 He are produced. The weak reaction rates are computed by integrating neu- 
trino distribution functions which are obtained by solving Boltzmann equations numerically. 
Using the neutrino distribution functions, the six weak interaction rates between neutron 
and proton are represented by 



rQ-m e 

Jo 

/•CO 

T ne +^ pl7e = K dp Ve 

JQ+m e 



ip VB - q) 2 - m\{Q - P^ 1+e ;::_ QyTy (i - um) 



(p. - QY - ml( P , e - Q) eiP ^ + 1 (i - UM) 



(25) 
(26) 



the reheating temperature decreases. This is because neutrinos are produced by the annihilation 
of electrons-positron pairs whose mean energy per one particle is larger than that of neutrinos. 
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K / dp v 
Jo 



(Pu e +Q) 2 -ml(p Ue + Q) 
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I _|_ e -(p, e +Q)/T 7 
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dp v 



(Pve ~ Q) 2 ~m 2 e (Q -Pu e ) 



UXPv.) 



pe~—>nv e 



roo 

K / dp v 
Jo 



(p„ e + Q) 2 - ml(Q + p Ve 



e ~{p^-Q)l T -i 4- 1 

2 

1 - fvAPve)) 



PI 



e (Pv e +Q)/ T l -f l 



r 



pv e ^rne 



K 



Q+m e 



dp v 



(Pve ~ QY ~ m l(Q - Pv 



1 4- e ~{P Ve ~Q)l T -i 



fvAPv 



(27) 
(28) 
(29) 
(30) 



where Q = m n — m p = 1.29 MeV and K is a normalization factor which is determined by the 



neutron life time r n as K ~ (1.636r n ) _1 and r n is obtained by the experiments [|27|]. From 
the above equations we can see that if neutrino and anti-neutrino distribution functions 
are decreased, both (3 decay rates T n ^ p = T n ^ pe -^ e + Y ne +^ v ^ + F and T r 



nue^pe 



pe u e - 



pe 



+ r 



pv e — me 1 



are simultaneously decreased by the following reasons. The 



dominant effects by the deficit of the distribution functions are to decrease the rates T r 



->pe 



Tpe-y^n and T p „ e ^ ne + which have the neutrino or anti-neutrino in the initial state. On the 



other hand, though the other rates r n ^ pe - l7e , T ne +^ pi7e and T pe -^ nUe which have the neutrino 
or anti-neutrino in the final state are slightly increased due to Fermi-blocking factor (l — f v ), 
the ratio of the difference Af u to (1 — f v ) is much smaller than that of Af u to f v , i.e. 



\Af v /(l - f v )\ « \Af v /f v \ for f v « 1. 



(31) 



Therefore, the enhancement is small and the latter effect is not important. In total, both 
weak interaction rates T n ^, p and T p ^ n are decreased and become smaller than those of SBBN. 
In Fig. H the weak interaction rates T n ^ p and T p ^ n are plotted. The solid lines denote the 
case of Tr = 10 MeV which corresponds to the standard big bang scenario. The dotted lines 
denote the case of Tr = 1 MeV. In the plot we can see that the insufficient thermalization 
of the neutrino distributions derives the changes of the weak interaction rates. 

The decrease of weak interaction rates gives significant effects on the abundance of 4 He. 
When the weak interaction rate T n ^ p decreases, Hubble expansion rate becomes more rapid 
than that of the weak interaction rate earlier. Namely the freeze-out time becomes ear- 
lier. Then the freeze-out value of neutron to proton ratio becomes larger than in SBBN 
and it is expected that the predicted 4 He abundance becomes larger. The above effect is 
approximately estimated by 



ay ~ +o.i9(-Ar riMp /r, 



(32) 



In Fig. ^ we plot the time evolution of the neutron to proton ratio. In Fig. ||(a) we 
change only the number of neutrino species in SBBN. The dotted line denoted the curve 
of iV° ff = 1.37 which corresponds to the effective number of neutrino species in the case of 
Tr = 2 MeV in the late-time entropy production scenario. Then we find that the predicted 
n/p curve is lower than that of Nf 1 = 3 due to only the speed down effects or the later 
decoupling. In Fig. [j](b) we plot the time evolution of n/p when we change the reheating 
temperature in the late-time entropy production scenario. The dotted line denotes the case 
of Tr = 2 MeV. Comparing to the case of iV^ ff = 1.37 in Fig. |6](a), the n/p ratio becomes 
larger. It is because the weak interaction rates are decreased by the deficit of the distribution 
function. Moreover in the case of Tr = 1 MeV the n/p ratio becomes much larger. 
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C. Neutrino thermalization and light element abundances 



Next we perform Monte Carlo simulation and the maximum likelihood analysis |28] to 
discuss how the theoretical predictions with the low reheating temperature scenario agree 
with the observational light element abundances. 

In Fig. [7| we plotted 4 He mass fraction Y as a function of Tr at 77 = 5 X 10~ 10 (solid 
line). The dashed line denotes the virtual 4 He mass fraction computed by including only 
the speed down effect due to the change of the effective number of neutrino species which is 
shown in Fig. |j. The dotted line denotes the predicted value of Y in SBBN at rj = 5 x 10~ 10 . 
For Tr > 7MeV, the solid line and dashed line are quite equal to the value in SBBN. As 
Tr decreases, both the solid and dashed lines gradually decrease because of the speed down 
effect due to the change of Nf 1 . The dashed line continues to decrease as the reheating 
temperature decreases. 

On the other hand, for T^ <j 2 MeV the effect that the weak interaction rates are 
weakened due to the deficit of the neutrino distribution function begins to become important 
and the predicted value of Y begins to increase as T R decreases. For T R < 1 MeV, since 
it is too late to produce enough electrons whose mass is about m e = 0.511 MeV, the weak 
interaction rates are still more weakened and Y steeply increases as Tr decreases. 

In Fig. [S] we plot the contours of the confidence level in the r]- Tr plane. The solid line 
denotes 95 % C.L. and the dotted line denotes 68 % C.L. The filled square is the best 
fit point between the observation and theoretical predictions. The observational data are 
consistent with the high baryon to photon ratio, rj ~ (3 — 6) x 10~ 10 . From Fig. |8| we find that 
Tr ^ 0.7 MeV is excluded at 95 % C.L. In other wards Tr as low as 0.7MeV is consistent 
with BBN. Then can be as small as 0.1 and it definitely influences the formation of the 
large scale structure and CMB anisotropy as is seen in SecJV|. 



V. HADRON INJECTION BY MASSIVE PARTICLE DECAY 
A. Hadron Jets and e + e~ collider experiments 

In the previous section we discussed only the case in which the parent massive particle 
4> decays into photons or the other electro-magnetic particles. In this section we consider 
the entropy production process along with the hadron injection, i.e. the case in which the 
massive particle has some decay modes into quarks or gluons. Then the emitted quark- 
antiquark pairs or gluons immediately fragment into hadron jets and as a result a lot of 
mesons and baryons, e.g. pions, kaons, nucleons (protons and neutrons ) are emitted into 
the electro-magnetic thermal bath which is constituted by photon, electron and nucleons. 

For example, if the gravitino ip^ is the parent particle which produces the large entropy, 
it could have a hadronic decay mode {e.g. ip^ — > ^qq ) with the branching ratio Bh — 0(a) 
at least even if the main decay mode is only ip^ — > 77(7 : photino) |29|. Then about 0.6 - 3 
hadrons are produced for ~ 1 — 100 TeV. In addition the emitted high energy photons 
whose energy is about m^/2 scatter off the background photons and could also produce the 
quark- ant iquark pairs through the electromagnetic interaction. For the cosmic temperature 
~ 0( MeV), the energy in the center of mass frame is \/s ~ 2 — 20 GeV for ~ 1 — 100 TeV. 
Then the number of the produced hadrons is about 2-7 which effectively corresponds to 
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the hadron branching ratio Bh ~ 10~ 2 if we assume that the hadron fragmentation is similar 
to the results of e + e~ collider experiments. Thus Bh should not become less than about 
10~ 2 for gravitino decay. [] For the other candidate, if the "flaton" is the parent particle as 
in thermal inflation model, it would also have a hadronic decay mode (0 — > gg ) || if the 
flaton mass is larger than 1 GeV. 

If once such hadrons are emitted to the electro-magnetic thermal bath in the beginning 
of BBN epoch (at T ~ lOMeV — 0.1 MeV), they quickly transfer all the kinetic energy 
into the thermal bath through the electro-magnetic interaction or the strong interaction. 
Through such thermalization processes the emitted high energy hadrons scatter off the 
background particles, and then they induce some effects on BBN. Especially, the emitted 
hadrons extraordinarily inter-convert the ambient protons and neutrons each other through 
the strong interaction even after the freeze-out time of the neutron to proton ratio n/p. For 
the relatively short lifetime (r^ ~ 10~ 2 sec — 10 2 sec) in which we are interested, the above 
effect induces the significant change in the previous discussion. Namely protons which are 
more abundant than neutrons are changed into neutrons by the hadron-proton collisions 
and the ratio n/p increases extremely. Because 4 He is the most sensitive to the freeze out 
value of n/p, the late-time hadron injection scenario tends to increase Y p . 

Reno and Seckel JJIJ investigated the influences of the hadron injection on the early 
stage of BBN. They constrained the lifetime of the parent particle and the number density 
comparing the theoretical prediction of the light element abundances with the observational 
data. Here we basically follow their treatment and apply it to the scenario of late-time 
entropy production with hadron injections. 

The emitted hadrons do not scatter off the background nucleons directly. At first hadrons 
scatter off the background photons and electrons because they are much more abundant than 
nucleons. For t < 200 sec, the emitted high energy hadrons are immediately thermalized 
through the electro-magnetic scattering and they reach to the kinetic equilibrium before they 
interact with the ambient protons and neutrons. Then we use the threshold cross section 
(av)^ ! i _ iN , for the strong interaction process N + Hi — > N' + ■ ■ ■ between hadron Hi and the 
ambient nucleon N, where N denotes proton p or neutron n. The strong interaction rate is 
estimated by 

r^Uiv = n N (av)^ N , 

where n^ is the number density of the nucleon species N, n is the baryon to photon ratio 
(= ns/ny), n B denotes the baryon number density (= n p +n n ) and Jn is the nucleon fraction 
(= n^/ns)- This equation shows that every hadron whose lifetime is longer than O(10~ s ) sec 
contributes to the inter-converting interaction between neutron and proton at the beginning 
of BBN. Hereafter we will consider only the following long-lived hadrons, (mesons, n , 
and Kl, and baryons, p, p, n, and n). For the relevant process ( N + — ► N' ■ --, and 



6 If the decay mode ipn — > gg (g : gluino) is kinematically allowed, the hadronic branching ratio 
becomes close to one. 
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iV + K~ — > AT' • ••, etc.), we can obtain the cross sections in pO| , pT| . Here we ignore K + 
interaction because n + K + — > p + if is the endothermic reaction which has Q = 2.8 MeV. 
We estimate the average number of emitted hadron species Hi per one (ft decay as 

N* = B h N jet f Hi ^, (34) 

where (N c h) is the averaged charged-particle multiplicity which represents the total number 
of the charged particles emitted per two hadron jets, /# 4 is the number fraction of the hadron 
species Hi to all the emitted charged particles, is the branching ratio of the hadronic 
decay mode and Nj et is the number of the produced jets per one (f> decay. 

Here it is reasonable to assume that the averaged charged particle multiplicity (N c h) is 
independent of the the source because the physical mechanism which governs the production 
of hadron jets is quite similar and does not depend on the detail of the origin only if the high 
energy quark-antiquark pairs or gluons are emitted. We adopt the data which are obtained 
by the e + e~ collider experiments. LEPII experiments (ALEPH, DELPHI, L3 and OPAL) 
recently give us the useful data for y/s = 130 — 172 GeV [pTp . We adopt the following fitting 
function for -y/s = 1.4 - 172 GeV 



(N ch ) = 1.73 + 0.268 exp ^1.42^/m(s/A 2 )J , (35) 

where y/s denotes the center of mass energy, the functional shape is motivated by the next-to- 
leading order perturbative QCD calculations, A is the cut-off parameter in the perturbative 
calculations and we take A = 1 GeV. In Fig. |9] we plot the charged particle multiplicity 
for y/s = 1 GeV — 100 TeV. The error of the fitting is about 10%. Using the available 
data [27,31,32], we obtain the emitted hadron fraction = n Hi / (A^), 



f w+ = 0.64, U- = 0.64, 
f K+ = 0.076, f K - = 0.076, 7^ = 0.054 (36) 
f P = fp = 0-035, f n = h= 0.034, 

where n Hl is the number of the emitted hadron species Hi which is defined as the value after 
both Kg and A had completely finished to decay. [] As we find easily, almost all the emitted 
particles are pions which are the lightest mesons. To apply the data of the e + e~ collider 
experiments, we take y/s = 2Ej et in (A^) where Ej et denotes the energy of one hadron jet 
because the (N c h) is obtained by the result for two hadron jets. 



Although the summation of is obviously more than one, it is because the experimental fitting 
of {N c h} is defined as a value before K$ and A begin to decay [^2|. Here we assume that fu. do 
not change significantly in the energy range y/s ~ 10 GeV - 100 TeV. Since we do not have any 
experimental data for the high energy region more than about 200 GeV, we extrapolate (Nch) to 
the higher energy regions and we take as a constant. 
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B. Formulation in Hadron Injection Scenario 



In this section we formulate the time evolution equations in the late-time hadron injection 
scenario. As we mentioned in the previous section, the hadron injection at the beginning 
of BBN enhances the inter-converting interactions between neutron and proton equally and 
the freeze out value of n/p can be extremely increased. Then the time evolution equations 
for the number density of a nucleon N(= p, n) is represented by 



dn 



N 



dt 



+ 3H(t)n N 



dn 



N 



dt 



N — >N 



,-K 



(37) 



weak 



where H(t) is Hubble expansion rate, [dn^ / dt] wea k denotes the contribution from the weak 
interaction rates which are obtained by integrating the neutrino distribution functions as 
discussed in Sec. [IV], see Eqs.fpB] - ^0|) , = p^/m^ is the number density of 
denotes the average number of the transition N —>■ N' per one <fr decay. 
The average number of the transition iV — > N' is expressed by 
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K 
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(38) 



where Hi runs the hadron species which are relevant to the nucleon inter-converting reac- 
tions, N Hi denotes the average number of the emitted hadron species Hi per one decay 
which is given by Eq. fl34j) and R^_^, N , denotes the probability that a hadron species Hi 
induces the nucleon transition N — ► N', 
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t h ^ is the decay rate of the hadron Hi and T^ s 
N , is the total absorption rate of Hi. 
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C. Hadron injection and BBN 

In this subsection we compare the theoretical prediction of the light element abundances 
in the hadron injection scenario to the observational light element abundances. In the 
computations we assume that the massive particle decays into three bodies (Ej et = m^/3) 
and two jets are produced at the parton level (Nj et = 2) In the computing we take the 
branching ratio of the hadronic decay mode Bh = O(10~ 2 — 1). 



The above choice of a set of model parameters Ej e t and Nj e t is not unique in general and is 
obviously model dependent. However, since {N^) has the logarithmic dependence of Ej e t, we 
should not be worried about the modification of Ej e t by just a factor of two so seriously. On 
the other hand in Eq. (^), the second term in the right hand side scales as cx Nj^/m^. For 
the modification of Nj et , therefore, we only translate the obtained results according to the above 
scaling rule and push the responsibility off onto 
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As we noted in the previous subsections, it is a remarkable feature that the predicted 
Y p tends to increase in the hadron injection scenario because 4 He is the most sensitive to 
the freeze-out value of the neutron to proton ratio. Since protons which are more abundant 
than neutrons are changed into neutrons through the strong interactions rapidly, the freeze 
out value of n/p increase extremely if once the net hadrons are emitted. In Fig. [10] we plot 
the predicted 4 He mass fraction 7 P as a function of Tr for (a) m^lOO TeV and (b) 
= 10 GeV. The solid curve denotes the predicted Y p . Here we take the branching ratio of 
the hadronic decay mode as B h = 1 (right one) and B h = 0.01 (left one). The dot-dashed 
line denotes B h = 0. The dashed line denotes the virtual value of Y p computed by including 
only the speed down effect due to the change of the effective number of neutrino species. 
The dotted line denotes the prediction in SBBN. 

As we mentioned in the previous section, the speed down effect due to deficit of the 
electron neutrino distribution function are not important for Tr > 7MeV. In addition since 
it is high enough to keep n/p ~ 1 for the cosmic temperature T > 7MeV, the enhancements 
of the inter-converting interaction between n and p due to the hadron emission do not 
induce any changes on the freeze-out value of n/p. As Tr decreases (Tr < 7MeV), Y p also 
decreases gradually because the speed down effect on the freeze-out value of n/p begins to be 
important. On the other hand, if a lot of hadrons are emitted when the cosmic temperature 
is T < 6 - 7 MeV and the ratio n/p is less than one, they enhance the inter-converting 
interactions more rapidly. As a result, the ratio n/p attempts to get closer to one again 
although the cosmic temperature is still low. Thus the above effects extremely increase the 
freeze-out value of n/p and is much more effective than the speed down effects. Namely the 
produced Y p becomes larger very sensitively only if Tr is just a little lower than 6-7 MeV. 
One can obviously find that this effect becomes more remarkable for the larger B^. 

To understand how it depends on mass, it is convenient to introduce the yield variable 
Y$ which is defined by 

= n^/s, (40) 

where s denotes the entropy density in the universe. Because Y^ is a constant only while 
the universe expands without any entropy production, it represents the net number density 
of 4> per comoving volume. For the simplicity let's consider the instantaneous decay of <j) 
and assume that the reheating process has been completed quickly. Because the radiation 
energy in the thermal bath or entropy s = 2iT 2 g*/4:5T R ' is produced only from the decay 
products of (f>, Y^ is approximately estimated using Tr and by 

Ys ~ 0.28—. (41) 

m<f, 

From the above equation, we can see that for the fixed value of Tr the net number of 0, i.e. 
the net number of the emitted hadrons, becomes larger for the smaller mass. Comparing 
Fig. [l(] (a) with Fig. (b), we find that the theoretical curve of Y p for the case of = 
10 GeV is enhanced more steeply and the starting point to increase Y p becomes higher than 
for the case of = 100 TeV. 

Since the other elements (D and 7 Li) are not so sensitive as 4 He, it is expected that 
the observational value of Y p constrains Tr most strongly. In order to discuss how a low 
reheating temperature is allowed by comparing the theoretical predictions with observational 
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values (D, 4 He and 7 Li), we perform the Monte Carlo simulation and maximum likelihood 
analysis as discussed in Sec. [IV]. In addition to the case of Sec. [IV] we take account of the 
following uncertainties, the error for the fitting of (N c h) as 10% |5T| and the experimental 
error for each cross section of the hadron interaction as 50%. Because there are not any 
adequate experimental data for the uncertainties of cross sections f30| , |3~l"]j , here we take the 
larger values to get a conservative constraint. 

In Fig. |TT] we plot the contours of the confidence level for = 100 TeV in the (t)-Tr) 
plane for (a) B h — 1 and (b) B h = 10 -2 . The solid line denotes 95 % C.L., the dotted line 
denotes 68 % CL. and the filled square is the best fit point between the observation and 
theoretical prediction for D, 4 He and 7 Li. The baryon to photon ratio which is consistent 
with the observational data is restricted in the narrow region, rj ~ (4 — 6) x 10 -10 . From 
Fig. [Il](a), we find that T R < 3.7 MeV is excluded at 95 % C.L. for B h = 1. On the other 
hand, from Fig. 0(b) we obtain the milder constraint that Tr < 2.5 MeV is excluded at 95 
% C.L. for Bh = 10~ 2 . In Fig. we plot the contours of the confidence level for = 10 
GeV in the same way as Fig. [11]. Compared to Fig. |TT , as we mentioned above, we find that 
the lower bound on the reheating temperature becomes higher for a smaller mass. From 
Fig. |12| we get the lower bound on the reheating temperature that Tr > 5.0 MeV (4.0 MeV) 
at 95 % C.L. for B h = 1 (B h = 10~ 2 ) 

In Fig. ^ the lower bound on T R as a function of are plotted for (a) B h = 1 and (b) 
B h = 10 -2 . The solid line denotes 95 % C.L. and the dotted line denotes 68 % C.L. As it is 
expected, the curve of the lower bound on Tr is a gentle monotonic decreasing function of 
m^. In Fig. |l3Ka), we can see that Tr should be higher than 4 MeV at 95 % C.L. for Bh = 
1 in = 10 GeV - 10 2 TeV0. On the other hand, in Fig. [15(b) we find that the constraint 
gets milder for B h = 10~ 2 . It is shown that T R < 2.5 MeV is excluded at 95 % C.L. for 
B h = 10~ 2 . In Fig. [| we find that N* s can be allowed as small as 2.8 for Bh = 1 ( 1.9 for 
B h = W- 2 ). 



D. Summary of hadron injection 

In this section we have seen that the BBN constraint on the reheating temperature 
becomes much more stringent if a massive particle has a branching to hadrons. For successful 
BBN the reheating temperature should be higher than 2.5 — 4 MeV for the branching ratio 
Bh = 1 — 10~ 2 . The hadron injection generally occurs if the late-time reheating is caused 
by the heavy particle with mass larger than ~ 1 GeV. Many candidates for the late-time 
reheating such as SUSY particles and flatons have such large masses and hence the constraint 
obtained here is crucial in constructing particle physics models based on SUSY or thermal 
inflation models. 

For the lower limit of the reheating temperature, the effective number of the neutrino 



9 Though we have adopted the experimental error of the each hadron interaction cross section as 
50% in the Monte Carlo simulation because of no data, the lower bound on Tr might become about 
10% higher than the above values if we adopt the more sever experimental error as 10% instead of 
50%. 
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species N^ s is given by 2.8 and 1.9 for E>h = 1 and 10~ 2 , respectively. Since the limiting 
temperature is close to the neutrino decoupling temperature, the deviation of N„ from the 
standard value (i.e. 3) is small and hence the detection may not be easy. 

However, from more general point of view, it is possible that light particles with mass 
1 GeV are responsible for the late-time reheating. In this case, as seen in the previous 
section, the reheating temperatures as low as ~ 0.7 MeV are allowed. For such low reheating 
temperature, neutrinos cannot be produced sufficiently Thus the effective number of the 
neutrino species N„ becomes much less than 3. This leads to very interesting effects on 
the formation of large scale structures and CMB anisotropies, which we discuss in the next 
section. 



VI. CONSTRAINTS FROM LARGE SCALE STRUCTURE AND CMB 

ANISOTROPY 

In this section, we discuss possibility to set constraints on the late-time entropy produc- 
tion from the large scale structure and CMB anisotropies. Hereafter, we only consider flat 
universe models with cosmological constant which are suggested by recent distant Super- 
novae (SNe) surveys [j33| , j34|| and measurements of CMB anisotropies . 



The late-time entropy production influences formation of the large scale structure and 
CMB anisotropies since the matter-radiation equality epoch is shifted if the effective number 
of neutrino species changes. The ratio of neutrino density to black-body photon density is 
Pvl Pj = (7/8)(4/ll) 4//3 ^. Therefore the redshift of matter-radiation equality can be written 
as a function of N„ 



• V ■ 



( 7 / 4 \ 4 /3 \ - L 
l fz eq = 4.03xl0 4 fi /i 2 (l + g (jy) AU , (42) 



where Qo is the density parameter and h is the non-dimensional Hubble constant normalized 
by lOOkm/s/Mpc. 

Let us now discuss distribution of galaxies on large scales. For a quantitative analysis, 
we define the matter power spectrum in Fourier space as P{k) =< \5k\ 2 >, where 8k is the 
Fourier transform of density fluctuations and <> denotes the ensemble average. Hereafter, 
we assume the Harrison-Zel'dovich power spectrum, which is motivated by the inflation 
scenario, as an initial shape of the power spectrum, i.e., P{k) oc k. As fluctuations evolve in 
the expanding universe, the shape of the power spectrum is changed. One often introduces 
the transfer function T(k) to describe this modification of the initial power spectrum as 
P(k) = AkT(k) 2 , where A is an arbitrary constant. In case of standard cold dark matter 
(CDM) dominated models, Bardeen et al. [^6| found a fitting formula: 



T(k) = ln(1+ ^ 34g) [l + 3.89g + (16.1g) 2 + (5.46g) 3 + (G^lg) 4 ]" 1 / 4 , (43) 

where q = fc/f2o/i 2 Mpc _1 when the baryon density is negligible small compared to the total 
density. It is easy to explain why q is parameterized by Q h 2 . This is because CDM 
density fluctuations cannot evolve and stagnate during a radiation dominated era. Only 
after the matter-radiation equality epoch, fluctuations can evolve. Therefore the CDM power 
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spectrum has a peak which corresponds to the horizon scale of the matter-radiation equality 
epoch. In fact, the wave number of the horizon scale at the equality epoch can be written 



as fc eq = J 2O (l + z eq)Ho, where H is the Hubble constant at present, that is proportional 
to Qoh 2 . In the actual observations, distances in between galaxies are measured in the units 
of h~ l Mpc. Therefore to fit the observational data by the CDM type power spectrum, we 
usually introduce so-called shape parameter T s = Q h. It is known that we can fit the galaxy 
distribution if r s ~ 0.25 ± 0.05 [[37] which suggests a low density universe. If the late-time 
entropy production takes place, however, we need to take into account dependence of 
the matter-radiation equality epoch (Eq. fl4"2"|)). Therefore T s should be written as 



r s = imn h/(i + o.227Nf). (44) 

We plot contours of T s on Qq — N„ plane in Fig. H. It is shown that smaller f2 is preferable 
for N„ < 3 with the same value of r s . We also plot the power spectra for Q = 0.3 and 
h = 0.7 with different iV^ ff 's in Fig. |15|. Here we do not simply employ the fitting formula 
but numerically solve the evolution of density fluctuations ||38|| . It is shown that the peak 
location of a model with smaller jV^ ff shifts to the smaller scale (larger in k) since smaller 
N„ makes the equality epoch earlier which means the horizon scale at the equality epoch 
becomes smaller. We have hope that current larege scale structrue surveys such as 2DF and 
Sloan Digital Sky Survey (SDSS) may determine the precise value of F s . 

Besides the shape of the power spectrum, the amplitude is another important observa- 
tional quantities to test models. On very large scales, the amplitude of the power spec- 
trum is determined by CMB anisotropies which are measured by COBE/DMR ||43|| . Since 
COBE/DMR scales are much larger than the horizon scale of the matter-radiation equality 
epoch, however, it is not sensitive to the transfer function T(k) but the overall amplitude A. 
In order to compare the expected amplitude of the power spectrum from each CDM model 
with large scale structure observations, we employ the specific mass fluctuations within a 
sphere of a radius of 8/i~ 1 Mpc, i.e., a 8 which is defined as 

al = < {5M/M{R)) 2 > R=8h -i Mp c 

= 2^2/ dkk 2 P(k)W(kR) 2 \ R=8h -i Mpc , (45) 

where W(kR) is a window function for which we employ a top hat shape as W(kR) = 
3 (sm(kR) — kR cos(kR)) /(kR) 3 . Eke et al. |39] obtained the observational value of a 8 



which is deduced from the rich cluster abundance at present as 

a 8 = (0.52 ± 0.0A)n o °- 52+0A3no . (46) 

Other estimates of erg 0] are agreed with their result. For CDM models with standard 
thermal history, the value of a 8 is a function of Qq an d h. With the late-time reheating, 
however, <jg for fixed f2o and h becomes larger. The reason is following. Since we fix Qq 
and h, the normalization factor A is same regardless of the value of N„. As is shown in 



Fig. 15, the amplitude of the power spectrum on 8/i -1 Mpc, i.e., as, is larger for smaller 
Nl E . In Fig. [TB|, we show the allowed region on the Q — h plane for Nff = 0.5, 2 and 3 for 
COBE-normalized flat CDM models with the Harrison-Zel'dovich spectrum. The shaded 



region satisfies the matching condition with the cluster abundance Eq. (|4q). For fixed h, 
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models with smaller N* s prefer lower Qq. Recently, the HST key project on the Extragalactic 



Distance Scale has reported that h = 0.71 ±0.06 (la) by using various distant indicators [41 



From SNe measurements, fio — 0.28 ± 0.8 for flat models (see Fig. 7 of Plfl). CDM models 
with N„ = 0.5 ~ 3 are still consistent with above value of h and Qq. However we expect 
further precise determination of Qo, h (from distant SNe surveys and measurements of CMB 
anisotropies) and ag (from 2DF or SDSS) will set a stringent constraint on N„ and T R in 
near future. 

Finally we discuss the CMB constraint on T R . Let us introduce temperature angular 
power spectrum Ci where I is the multipole number of the spherical harmonic decomposition. 
The rms temperature anisotropy of CMB can be written as < \ AT/T\ 2 >= X^(2£+l)C£/47r. 
Using Ce, we can extract various important information of cosmology, such as the curvature 
of the universe, flo, cosmological constant, h and so on (see, e.g., J42|). In fact, we can 
measure the matter radiation equality epoch by using the height of peaks of Cg. The peaks 
are boosted during the matter-radiation equality epoch. If the matter-radiation equality is 
earlier, the correspondent horizon scale is smaller. Therefore we expect lower heights for first 
one or two peaks since these peaks are larger than the horizon scale at the equality epoch 



and do not suffer the boost as is shown in Fig. [L7|. With the present angular resolutions 
and sensitivities of COBE observation [13] or current balloon and ground base experiments, 
however, it is impossible to set a constraint on N„. It is expected that future satellite ex- 
periments such as MAP || and PLANCK [Hj will give us a useful information about N„. 
From Lopez et al.'s analysis 0, MAP and PLANCK have sensitivities that 5N° S > 0.1 
(MAP) and 0.03 (PLANCK) including polarization data, even if all cosmological param- 
eters are determined simultaneously (see also Fig. [T7|). From such future observations of 
anisotropies of CMB, it is expected that we can precisely determine T R . 



VII. CONCLUSION 

In this paper we have investigated the various cosmological effects induced by the late- 
time entropy production due to the massive particle decay. The neutrino distribution func- 
tions have been obtained by solving the Boltzmann equations numerically. We have found 
that if the large entropy are produced at about t ~ 1 sec, the neutrinos are not thermalized 
very well and hence do not have the perfect Fermi-Dirac distribution. The deficits of the 
neutrino distribution functions due to the insufficient thermalization decrease the Hubble 
expansion rate and weakens the weak interaction rates between proton and neutron. The 
above two effects changes the freeze-out value of n/p significantly. Especially the produced 
4 He mass fraction Y is so sensitive to n/p that the predicted value of Y is changed drastically. 
Comparing the theoretical predictions of D, 4 He and 7 Li to the observational data, we have 
estimated the lower bound on the reheating temperature T R after the entropy production. 
We have found that T R < 0.7 MeV is excluded at 95 % C.L. In other wards, T R can be as 
low as 0.7 MeV. Then the effective number of neutrino species Nf 1 can be as small as 0.1. It 
is enough sensitive for the ongoing large scale structure observations such as 2DF and SDSS 
or future satellite experiments (MAP and PLANCK) of CMB anisotropies to detect such 
modifications on N„ and we can find out the vestige of the late-time entropy production. 

Furthermore, we have also studied the case in which the massive particle has some decay 
modes into quarks or gluons. In this scenario, a lot of hadrons, e.g. pions, kaons, protons 
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and neutrons, which are originated by the fragmentation of the high energy quarks and 
gluons are injected into thermal bath. The emitted hadrons extraordinarily inter-convert 
the ambient protons and neutrons each other through the strong interaction even after the 
freeze-out time of the neutron to proton ratio n/p. Then the predicted value of Y increase 
extremely and we can constrain Tr and the branching ratio of the hadronic decay mode 
Bh comparing to the observational light element abundances. We have found Tr should be 
higher than 2.5 MeV - 4 MeV at 95 % C.L. for B h = KT 2 - 1. The above results tells us 
that N„ can be as small as 1.9 - 2.8 even in the hadron injection scenario for = 10~ 2 - 
1. Then it still may be possible to detect the modifications on iV^ ff by MAP and PLANCK. 
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APPENDIX: REDUCTION OF COLLISION INTEGRAL 



This appendix shows how we can reduce the nine dimensional integrals in Eq. (^) of 
the collision term Cj iScat for the scattering process into one dimensional integrals. Notice 
that, since we treat the massless neutrino, the norm of the neutrino momentum equals to 
its energy \pi\ = E^. Here we divide the collision term into two parts: 

C i , acat = -F + B 1 (Al) 

where F represents the forward process and B represents the backward process. They are 
given by 

F =wJ isS? / I + — 3 - ^ |M|2A - (A2) 



B =^- [ d?2 - [ dPs - [ dPA ,A 2nY5\ Pl + P 2-p 3 - Pi )S\M\ 2 A B , (A3) 

2E 1 J 2£ 2 (2tt) 3 J 2E 3 (2ny J 2£ 4 (27r) 3V ' V 

where g e = 2 and the phase space factors are given by 

A F = /i(£i)/ 2 (£ 2 ) (1 - fc(E 3 )) (1 - U(E A )) , (A4) 

A B = (1 - (1 - / 2 (£ 2 )) fi(E 3 )f 2 (E 4 ). (A5) 

The integral over d 3 p4 is immediately done using 5 3 (pi + P2 — Ps — Pa)- From the 
momentum conservation, |p 4 | is given by 

|p 4 | 2 = E\ = El + 2E 2 Rcosr] + R 2 , (A6) 
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where R = p^ — p 3) R = \R\ and cost] = R • p 2 /(|P2|R)- 

The remaining delta function 5(Ei + E 2 — E 3 — E 4 ) shows the energy conservation low 
which is given by 



E\ = E\ + E\ + E\ + 2{E\E 2 — E\E 3 — E 2 E 3 ) 
We can generally take the momentum axes as 

R= (0,0, R) 

p 2 = (E 2 sin 77 sin <fi, E 2 sin 77 cos 0, E 2 cos 77) , 
p 3 = (E 3 sin f, 0,£ 3 sin 0, 

where 

77i2 tti2 p2 



cos£ 



Then | cos£| < 1 demands 



\E\ — E 3 \ < R < i?i + £^3. 



The volume element of p 2 is given by c?p| = E 2 dcosr]d(j) and from Eqs. 
azimuthal angle is obtained by 



cos 77 = 

Then | cos 77 1 < 1 demands 



R 2 — (Ei — E 3 ) 2 — 2E 2 (Ei — E 3 
2E 2 R 



\E\ — E 3 \ < i? < + 2i?2 — -^3- 
From Eqs. ( A12 ) and ( A14 ), we can obtain the allowed region of R, 
\Ei - E 3 \<R< Inf[£i + E 3 , E x + 2E 2 - E 3 \. 



(A7) 

(A8) 
(A9) 
(A10) 

(All) 

(A12) 
) and ( |A7|) the 

(A13) 
(A14) 



(A15) 



Since the volume element of p 3 is given by dpi 



2irE 2 dE 3 d cos 9 where cos 9 



Pi ' P3/ '(E 3 Ei), the differential angle element is evaluated by 

R 



dcos9 



2E X E 3 



-dR. 



(A16) 



From Eq. ( |A15| ) we can see that the integration can be performed in the four allowed 
intervals, 



F + B 



128 
1 



— ^ / dE 3 / dE 2 / dR / /|M| 2 (-A F + A B ) 
iEf Jo Jo J Jo 2n 



128£f 



El 



+ 



X 



2tt 



dE, 
dE 3 



E 3 



dE, 



E 1 +2E 2 -E 3 



dR + 



E 3 



dE 9 



E1-E3 

Ei+2E 2 -E 3 



Ei 



dE, 



E a 



Ei+E 3 



dR 



E1+E3 



E1+E3 



dR + 



dE 3 



Ei 



dE 9 



Ei-Ez 

Ei+E-i 



E 3 



dR 



E1+E3 



f n ^S\M\\-A F + A B ) 

JO Z7T 



(A17) 
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Even though we only show the case of of u e here, we can get the same procedure for and 
v T if Cy and Ca are replaced by Cy and Ca- As we also noted in Sec fj ; we assume that 
electrons obey the Boltzmann distribution function e~ E l T . In addition, since neutrinos are 
massless, the energy momentum conservation gives pi ■ p^ = p 2 ■ P3 in the elastic scattering 
process. The above assumptions simplify the integrations still more. 

For the forward reaction, v{p\) + e ± (p 2 ) —> ^(^3) + e ± (p4), the phase space factor is given 

by 



A F = f u (E 1 )(l-f u (E 3 ))exp[- 



T 



(A18) 



Then F\ and F2 in Eq. flT2"| ) are analytically estimated as 
^1 = 



E 3 f E 1 +2E 2 -E 3 roc rE 1 +E 3 

dE 2 / dR+ dE 2 dR 

JEi-E-i JE 3 JE 1 -E 3 



2tt 



S\M\ 2 e 



2tt 256(C 2 + C 2 A )G 2 F 



2T 4 



El + El + 2T(E 1 - E 3 ) + AT 4 



rji'l 



EfE 2 3 +2E 1 E 3 (E 1 + E 3 )T 



+2{E 1 + E 3 ) 2 T 2 + A(E 1 + E 3 )T A + 8T 4 



_£3 

e t , 



(A19) 



Fo = 



E 3 rE 1 +2E 2 -E 3 roo rE!+E 3 

dE 2 / dR+ dE 2 dR 

-£i+£ 3 J-Ei+_B 3 JE 3 J-_Bi+B 3 



2- 



5|M| 2 e-^ 



2T 4 (£ 2 + £ 2 - 2T{E 1 -E 3 )+ 4T 4 )e 



4 El -E 3 



T 



2tt 256(C 2 + C\ )G\ 
E 2 1 E 2 3 + 2E 1 E 3 (E 1 + E 3 )T 



+2(E 1 + E 3 ) 2 T 2 + A(E 1 + E 3 )T 6 + 8T 4 



e t . 



(A20) 



On the other hand, for the backward reaction, v(p\) + e ± {p 2 ) <— ^(^3) + e ± (p i ), the phase 
space factor is given by, 



As = (1 — f u {Ei)) f v {E 3 ) exp 1 



T 



(A21) 



Then we can analytically obtain Bi and B2 in Eq. dl2|) as 



E 3 r E 1 +2E 2 -E 3 /-oo 

d£ 2 / dR+ dE 2 I dR 

(I JEx-E 3 JE 3 JE 1 -E 3 



2-d0 g|M| 2 e~ t 3 



2 ; 



2T 4 (£ 2 + £ 2 + 2T(E 1 - £ 3 ) + 4T 4 )e" 



■Bl-g 3 



2tt 256(C£ + C^G 2 , 
E 2 E^ + 2E\E 3 {E\ + E 3 )T 



+2(E 1 + E 3 ) 2 T 2 + 4(Ei + £3)^ + 8T 4 



e T 



(A22) 



Bo 



E 3 rE 1 +2E 2 -E 3 poo rEx+E 3 

dE 2 / dR+ dE 2 dR 

-E!+E 3 J-E!+E 3 JE 3 J~Ex+E 3 



2n 



2tt 256(C 2 + C\)G 2 F 
2T 4 (El + E 3 2 - 2T(E l - E 3 ) + 4T 4 ) - T 2 [£ 2 £ 3 2 + 2E 1 E 3 (E l + £ 3 )T 



+2(Ei + E 3 ) 2 T 2 + 4(Ei + E 3 )T J + 8T 4 



(A23) 
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TABLES 



Process 


S|M| 2 


v e + e" 
v e + e+ 
v e + v e 




v e + e" 
f e + e+ 
e+ + e" 


32G$ [(C v 
32G$ [(C v 
32G 2 F [(C v 


+ C A ) 2 ( Pl 

-c A ) 2 ( Pl 

+ C A ) 2 ( Pl 


■P2? + (C V 

■P2? + (C V 

■Pa) 2 + (C v 


- C A ) 2 ( Pl 
+ C A ) 2 ( Pl 

-c A ) 2 ( Pl 


■PA?] 
■PA?] 
■PS?] 



TABLE I. Matrix elements for electron neutrino interactions. Gf is the Fermi coupling con- 



stant. Here we take Cv = \ + 2 sin 2 9w, C A = \ and the weak mixing angle sin 2 8\v — 0.231. 



Process 


S|M| 2 


+ e" 


-> + e~ 


32G|n 


(CV 


+ Ca?(pi 


■P2? + (Cv 


-Ca?(pi 


■PA? 




Vn + e + 


^ + e + 


32G$ 


(Cv 


-Ca?(pi 


■P2? + (C V 


+ Ca?(pi 


■PA? 






-» e+ + e~ 


32G| 


\Cv 


+ Ca?(pi 


■PA? + (C V 


-Ca?(pi 


■PS? 





TABLE II. Matrix elements for muon neutrino or tau neutrino interactions. Gf is the Fermi 
coupling constant. Here we take Cy = Cy — 1 = — \ + 2 sin 2 6\v, C A = C A — 1 = — \ and the weak 
mixing angle sin 2 6y/ ~ 0.231. 
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FIG. 1. Time evolution of the cosmic temperature (a) for Tr = 10 MeV, and (b) for 
Tr = 2 MeV. The dashed line denotes the neutrino temperature which can be defined only when 
they are thermalized sufficiently and have the perfect Fermi-Dirac distribution. 
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FIG. 2. Time evolution of the fraction of the energy density of v e (solid curve) and (dashed 
curve) to that of the standard big bang scenario for (a) Tr = 10 MeV and (b)Tjj = 2 MeV. Since 
the interaction of u T is as same as the curve of also represents v T . 



28 



Oh 



10 1 i — 

10° p 

10" 1 r 

10" 2 r 

IO- 3 r 
!(a) 



1 1 1 I I I I 



V 10' Me% 



J I 



J I I I I I 



10 1 
10° 
10" 1 

lO- 2 



10 



-4 



.1 


i 

PA 


1 


I 1 1 


1 1 MM 


T R = fe'jteVI 


i iiiiiiij i iiniiij 1 1 nun 1 1 
i 


1 1 MM 


' — — • . 

S \ — 
N \. - 
\ \. 

i i i i i NiV 


,i 


1 


l 



PA 



FIG. 3. Distribution function of v e (solid curve) and Vn (dashed curve) I 
and (b)for Tr = 2 MeV. The dotted curve is the Fermi-Dirac distribution 
interaction of v T is as same as v^,, the curve of also represents v T . 
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FIG. 4. Effective number of neutrino species N° s as a function of reheating temperature Tr. 
The top horizontal axis denotes the lifetime which corresponds to Tr. 



30 




10° 10" 1 lo- 2 
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FIG. 5. Weak interaction rates (sec -1 ) between neutron and proton. The upper curves are 
T n ^ p . The lower curves are r p _> n . The solid lines denote the case of Tr = 10 MeV which 
corresponds to the standard big bang scenario. The dotted lines denote the case of Tr = 1 MeV 
in the late-time entropy production scenario. Notice that T~\ p reaches r n = 887 sec in the low 
temperature. 
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FIG. 6. Evolution of neutron to proton ratio as a function of the temperature, (a)when we 
change only the number of neutrino species in the standard big bang scenario, and (b)when we 
change the reheating temperature in the late-time entropy production scenario. The dashed line 
is the thermal equilibrium curve (= e~®l T ). 
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FIG. 7. 4 He mass fraction Y p as a function of Tr (solid line) at r\ = 5 x 1CP 10 . The dashed line 
denotes the virtual 4 He mass fraction computed by including only the speed down effect due to the 
change of the effective number of neutrino species which is shown in Fig. ||. The dotted line denotes 
the value predicted in SBBN at r/ = 5 x 10~ 10 . The long-dashed line denotes the observational 2 
a upper bound, Y obs ~ 0.252 which is obtained by summing the errors in quadrature. The top 
horizontal axis represents the lifetime which corresponds to Tr. 
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FIG. 8. Contours of the confidence level in (77, Tr) plane. The inner (outer) curve is 68% (95%) 
C.L.. The filled square denotes the best fit point. The right vertical axis denotes the lifetime which 
corresponds to Tr. 
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FIG. 9. Plot of the charged particle multiplicity (N c h) for the center of mass energy 
lGeV- 100 TeV. 
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FIG. 10. Plot of the predicted 4 He mass fraction Y p as a function of Tr for (a) ni0=lOO TeV 
and (b) =10 GeV at r\ = 5 x 10~ 10 . The solid curve denotes the predicted Y p where we take 
the branching ratio of the hadronic decay mode as Bh = 1 (right one) and Bh = 0.01 (left one). 
The dot-dashed line denotes Bh = 0. The dashed line denotes the virtual Y p curve computed by 
including only the speed down effect due to the change of the effective number of neutrino species. 
The dotted line denotes Y p in SBBN. The long-dashed line denotes the rough observational two 
a upper bound that Y p should be less than about 0.252. The top horizontal axis represents the 
lifetime which corresponds to Tr. 
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FIG. 11. Contours of the confidence levels for = 100 TeV in (77, Tr) plane for the branching 
ratio of the hadronic decay mode (a) Bh = 1 and (b) Bh = 10 -2 . The solid line denotes 95 % 
C.L. and the dotted line denotes 68 % C.L. The filled square is the best fit point between the 
observation and theoretical prediction for D, 4 He and 7 Li. The right vertical axis represents the 
lifetime which corresponds to Tr. 



37 



10 



> 

CD 



:(a) 



i 1 — r 



II I I I 



::(b) 



m rf ,= 10GeV 



B h = 1 



1 

10 -io 



J 



i 1 — r 



FTT 



m =lOGeV 
B h = 10"* 



J I I 



CV2 
I 

o 



o 

CD 
GO 



- I 



O h 



10 



-9 



10 



-9 



7] 



7] 



FIG. 12. Contours of the confidence levels for = 10 GeV for the same theory parameters 



as in Fie. 11. 
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FIG. 13. Lower bound on Tr as a function of for the branching ratio of the hadronic decay 
mode (a) = 1 and (b) = 10~ 2 . The solid line denotes 95 % C.L. and the dotted line denotes 
68 % C.L. The right vertical axis represents the lifetime which corresponds to Tr. 
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FIG. 14. Contours of T s = 0.2 (bold), 0.3,0.4,0.5 and 0.6 on the (Sl ,Nf) plane for h = 0.7. 
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FIG. 15. Matter power spectra P(k) of CDM models with Nf = 0.5,2, and 3. We take 
Qo = 0.3, h = 0.7, and f^/i 2 = 0.02 where Qb is the baryon density parameter. 
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FIG. 16. Allowed region on — h plane from observational values of as deduced from the rich 
cluster abundance at present for flat CDM models. Models with N„ = 0.5, 2, and 3 are plotted. 
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FIG. 17. Power spectra of CMB anisotropics (left top panel) and polarization (right top panel) 
of models with =3,2 and 0.5. Bottom two panels show (Ci(N e ff) — C^(3))/C^(3) with 

N^ff = 2.9,2.5 and 2 for CMB anisotropies (left bottom) and polarization (right bottom). 
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